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ABSTRACT 


N. 


X.  v 

A  piecewise  linear  finite  element-based  method  of  lines  is  presented  for 
the  numerical  solution  of  coupled  parabolic  partial  differential  equations 
which  model  biological  and  physicochemical  reaction-diffusion  processes  in  one 
space  dimension.  The  vertical  lines  emanating  from  the  space  nodes  in  this 
method  change  at  automatically  selected  times  when,  in  order  to  control  a  norm 
of  the  space  discretization  error,  adaptive  spatial  regridding  occurs.  The 
regridding  algorithm  is  an  extension  of  one  described  previously  -by  tiye 
Authors  .[Tp  and  is  implemented  in  the  program  FEM0L1,  which  uses  the  LSODI 
package  of  lUndmarsh- and  Pointe^to  integrate  the  ordinary 

differential  equations  in  time  along  the  vertical  lines.  Computational 
results  show  that  the  method  is  efficient,  that  a  posteriori  estimates  of  the 
space  discretization  error  are  accurate,  and  that  the  adaptive  procedure 
reliably  controls  the  space  discretization  error. 
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1 .  INTRODUCTION 

Reaction-diffusion  processes  occur  in  many  branches  of  biology  and 
physical  chemistry.  Examples  include  substrate  transport  and  consumption  in 
the  microcirculation,  flame  propagation  in  combustion,  nerve  conduction,  and 
interactions  of  mobile  populations  in  ecosystems.  These  diverse  phenomena  are 
often  modeled  by  initial  boundary  value  problems  (IBVPs)  in  which  the 
governing  parabolic  partial  differential  equations  (PDEs)  are  nonlinear  and 
coupled  only  through  the  rates  at  which  physical  components  react.  Such  IBVPs 
in  one  space  dimension  are  the  problems  considered  in  this  paper. 

Solutions  of  these  problems  may  decay  to  steady  states,  oscillate  in  time 
or  evolve  as  localized  traveling  waveforms.  Spatial  rezoning  or  regridding  in 
numerical  methods  for  time-dependent  PDEs  with  the  latter  type  of  solutions 
has  become  quite  popular.  Grid  evolution  is  generally  governed  by  a  feedback 
procedure  in  the  "physical"  reference  frame  or  by  explicit  mappings  from 
"physical"  to  "computational"  coordinates,  either  in  an  a  priori  manner  or  as 
part  of  a  feedback  procedure.  Schemes  utilizing  coordinate  mappings  are 
surveyed  in  Thompson  et.  al.  [22]. 

In  this  paper  we  describe  the  implementation  of  a  variation  of  the 
classical  method  of  lines  (MOL)  in  which  the  space  grid  is  updated  in  a 
feedback  (adaptive)1  procedure.  The  MOL  reduces  an  IBVP  through  space 
discretization  into  an  initial  value  problem  for  a  system  of  ordinary 
differential  equations  (ODEs)  in  time.  The  ODEs  determine  spatial  parameters 
of  the  solution  on  vertical  lines  extending  in  time.  The  method  described 
here  is  based  on  a  finite  element  formulation  using  linear  elements  in 


*The  terms  "adaptive"  and  "feedback"  are  used  interchangeably  to  describe 
numerical  methods  here.  For  analyses  which  distinguish  between  the  two,  see 
Rheinboldt  [20]  and  Babuska  and  Vogelius  [1]. 


space.  We  refer  Co  this  approach  as  Che  FEMOL.  The  ODEs  in  Che  FEMOL 
determine  nodal  values  of  Che  solution.  They  are  integrated  via  variable- 
order,  variable-step  implicit  formulas  as  implemented  in  the  LSODI  package  of 
Hindmarsh  and  Painter  (cf,  [14],  [15]).  In  the  adaptive  FEMOL  grids  change 
discontinuously  in  time,  as  nodes  are  both  added  and  deleted.  Such  regridding 
is  often  said  to  be  locally  "static,”  in  contrast  to  "dynamic"  regridding, 
where  changes  occur  continuously  in  time. 

Most  methods  described  in  the  literature  which  employ  static,  dynamic  or 
combinations  of  both  types  of  regridding  can  be  viewed  as  MOL  extensions  or 
variations.  In  K.  Miller  and  R.  Miller  [16],  Gelinas  et .  al.  [12]  and  K. 
Miller  [17],  the  moving  finite  element  method  was  developed,  where  a  fixed 
number  of  nodes  are  used  and  nodal  positions  are  computed  together  with  the 
solution  values  at  the  nodes.  The  concept  of  moving  (in  time)  meshes  was  also 
used  in  Davis  and  Flaherty  [12]  and  Flaherty  et.  al.  [10].  Many  regridding 
methods  for  time-dependent  PDEs  have  attractive  properties.  The  reader  is 
especially  referred  to  Berger  and  Oliger  [4],  Dwyer  et .  al.  [9],  Gannon  [11] 
and  Harten  and  Hyman  [13]. 

While  these  methods  incorporate  various  principles  and  constraints  in 
feedback  approaches,  all  have  similar  objectives.  Regridding  is  expected  to 
yield  high  accuracy,  reliability  and  robustness  per  computational  cost.  It  is 
difficult  to  quantitatively  compare  various  feedback  (adaptive)  methods, 
despite  the  fact  their  goals  are  related,  but  it  should  be  possible  to  measure 
the  success  or  failure  of  each  in  achieving  its  goals.  Unfortunately,  the 
goals  often  are  not  well  formulated,  quantified  or  used  in  supporting 
computational  experiments. 

The  primary  goal  of  the  adaptive  FEMOL  is  to  control  the  space 
discretization  error  of  the  approximate  solution  as  measured  in  a  weighted 


L>2  gradient  norm  (weighted  norm) .  This  norm  arises  naturally  in 

connection  with  the  finite  element  formulation.  Error  control  is  achieved 
through  control  of  computed  error  estimates.  The  a  posteriori  estimator  used 
here  is  similar  to  that  analyzed  for  linear  parabolic  PDEs  by  the  authors  in 
[6],  [7].  The  estimator  is  obtained  by  summing  local  error  indicators.  These 
are  formed  from  PDE  residuals  evaluated  with  the  approximate  solution  inside 
each  of  the  elements . 

The  regridding  strategy  described  here  extends  an  earlier  version  of  the 
authors  [7]  and  was  summarized  for  linear  PDEs  in  Bieterman  [5].  A  grid  is 
retained  until  the  estimated  error  exceeds  a  preset  tolerance.  The  error  is 

lowered  below  this  tolerance  by  adding  and  deleting  nodes.  The  most  important 

part  of  grid  construction  is  a  pattern  recognition  procedure  used  to  determine 

the  grid's  "shape."  In  this  procedure  information  is  extracted  from  the  local 

error  indicators,  reduced  to  a  grid-independent  form,  and  relevant  features 
taken  from  this  reduced  data  are  compared. 

This  paper  is  organizied  as  follows.  Section  2  contains  mathematical 
details  of  the  problem  and  classical  version  of  the  method.  The  salient 
features  of  the  adaptive  FEMOL  are  described  in  Section  3,  along  with  the 
method's  goals  and  the  basic  strategy  used  to  achieve  them.  The  achievement 
of  the  goals  depends  first  and  foremost  on  the  quality  of  the  error 
estimator.  The  estimator  and  the  local  indicators  used  to  form  it  are  the 
subjects  of  Section  4.  In  Section  5  we  introduce  the  notions  of  "shape"  and 
"intensity"  to  describe  a  grid.  The  specific  strategy  used  to  construct  a 
grid  is  to  somewhat  directly  control  these  two  grid  properties.  This  strategy 
and  details  of  the  refinement/derefinement  algorithm  are  given  in  Section  5. 

In  Section  6  we  describe  the  selection  of  the  parameter  which  controls  grid 
intensity.  The  pattern  recognition  procedure  used  to  select  the  function 


controlling  grid  shape  is  presented  in  Section  7.  Section  8  contains  the 
results  of  many  carefully  conducted  computational  experiments  with  four 


I 


reaction-diffusion  problems.  These  results  enable  one  to  quantitatively 
evaluate  the  method's  performance.  They  show  how  the  regridding  strategy  is 
carried  out,  how  the  error  estimator  works,  and  the  dependence  of  error 
control  on  error  estimation.  The  final  section  summarizes  many  of  the 
important  aspects  of  the  method. 
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2.  PROBLEM  AND  CLASSICAL  METHOD  DESCRIPTION 

Let  fl  -  (x:  —  <  <  x  <  30+  <  •},  3Q  -  3fl~  U  3Q+  and  ii  ■  il  U  38. 

We  are  interested  in  finding  u  »  {u*( t ,x) } jNpde  ^or  c  ^  1 0 * T 1  and 
x  €  ft  which  satisfies  the  differential  equations 

u*  -  (ai(x)  u^)x  -  f*(t ,x,u) ;  x  €  a,  t  6  (0,T],  i  -  l.NPDE,  (l.a) 

the  boundary  conditions 

oj1  (x)  u1  +  0i(x)  ux  -  gi(t,x);  x  €  30,  t  €  (0,T],  i  =■  l.NPDE.  (l.b) 

and  the  initial  conditions 

u*  -  u*(x);  X€0,  t=0,  i  =  l.NPDE.  (l.c) 

For  x  f  0  it  is  assumed  that  a*(x)  >0,  i  »  l.NPDE,  with  a^Cx)  >  Sq  > 
0  for  at  least  one  index  i.  The  functions  a*  and  0*  determine  the  type 
of  boundary  conditions  (BCs).  It  is  assumed  that  0*'(3fi  )  <  0,  0*(3J2+)  >  0 
and  ct*(x)  >  0  for  x  €  3fl  and  i  *  l.NPDE.  With  reasonable  problem  data, 
the  solution  u  of  Eqs.  (1)  exists  and  is  a  smooth  mapping  of  [ 0 , T ]  into 
the  Hilbert  space 

H  -  {v  -  {v1}^!  >NpDE:  v1,^  €  L2(0);  i  -  l.NPDE}. 

Here,  l2(3)  denotes  the  usual  space  of  square  integrable  functions  on  0, 
with  norm  and  inner  product  '**q  an<*  <*,•>,  respectively.  Provided 

appropriate  conditions  on  {a*  ,a*  ,0^}^.j  NPDE  hold,  t*ie  semi-norm 

NPDE  .  .  ,  i/ 

I  I  Ml  I  5  {I  <a  vx*vi>}/2 

i-1 
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is  a  norm  on  H . 

We  will  denote  by 

6  -  {3Q-  “  xQ  <  <  x2  <  ...  <  xN_L  <  xN  -  3fl+}  (2) 

a  space  grid,  with  nodes  {xn}  and  elements  {(xn_^,xn)},  an<^  write 


h 


n 


xn  “  xn-l  for  n  -  l.N. 


(3) 


S(6)  c  H  denotes  the  finite  element  subspace  of  functions  which  are  linear  on 
each  element  of  6 . 

The  classical  version  of  the  FEMOL  is  formulated  in  the  present  setting  as 
follows.  With  given  5,  the  semi-discrete  approximation  U:  [0,T]  ♦  S(5) 
of  u  satisfies  the  equations 


<uJ(tf.)f«i>+<aiU^(t,.),«1>  +  ai(x)  4^  &i<t.x>*i<x>/!2+ 


8i(x) 


(4. a) 


<F1(t,.,l)),*i>  +  S-&L  gi(tfx)*i(x)/?jjl  ;  t  €  (0,T),  i  -  1, 


8  (x) 


NPDE, 


$  =.  {$J} 


j-l,NPDE 


€  S(6) 


and 


U(0,.)  -  0(0,.).  (4.b) 

U(0,«)  €  S(6)  is  the  linear  interpolate  of  Uq(»)  and  F(t,*,ft)  £  S(5) 
interpolates  f(t,*,U).  Here  it  has  been  assumed  for  simplicity  that  8i(x)  ^ 
0  for  all  i.  Otherwise,  S(6)  and  Eqs.  (4)  are  modified  by  usual  finite 
element  techniques. 

Eqs.  (4)  are  equivalent  to  requiring  that  the  vector  U[t]  - 
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{^(C  »xn))i,1  (ipug  satisfies  the  ODE 


i*l ,  NPDE 
n-0,N 


M  *  TT  U [ t]  +  A  •  U[t] 
at 


M  •  F[t,U]  +  B[t] ;  t  €  <0,T], 


(5.a) 


and  initial  condition 


U[0]  -  U[0]. 


(5.b) 


In  Eqs.  (5),  U [0 ]  and  F[t,U]  are  the  NPDE  •  (N+l )-dimensional  vec‘  s  of 

nodal  values  of  U(0,»)  and  F(t,«,U).  M  and  A  are  symmetric  matr  3 

which  are  positive  definite  and  nonnegative,  respectively.  With  an 
appropriate  ordering  of  the  nodal  values,  these  matrices  and  the  product  of 
M  and  the  Jacobian  of  F[t,U]  with  respect  to  U[»]  have  half  bandwidths 
2*NPDE  -  1  (i.e.  at  most  4*NPDE  -  1  nonzero  entries  per  row).  The  vector 

B[t]  has  at  most  2* NPDE  nonzero  entries  which  correspond  to  the  values 

(.1(*n)gl(t,*n)/81(«n))1,l>(IpDE  . 

n«0  and  N 

Recall  that  6*  ^  0  was  assumed.  If  6i(xQ)  *  0  and  c^Cxq)  f  0  (Dirichlet 
BC)  ,  for  example,  this  condition  is  implemented  by  modifying  entries  of  M, 

A  and  B[t]  to  obtain  the  equations 

°‘i(x0)  It  ^l(t‘xo)  "  gi(t,x0);  t  >  0,  and  Oi(0,x())  -Ui(0,xQ). 


The  FEMOL  approximate  solution  U(»,»)  is  obtained  by  numerically  solving 
Eqs.  (5).  This  is  accomplished  with  the  "stiff"  implicit  backward 
differentiation  formulas  in  the  LSODI  package  of  Hindmarsh  and  Painter  (cf. 
[14],  [15]).  LSODI  takes  an  input  time  discretization  error  tolerance  TOL  > 
0,  advances  with  internally  chosen  time  steps  and  integration  orders,  and 
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returns  U[t] ,  the  approximate  solution  of  Eqs.  (5)  and  the  vector  of  nodal 
values  of  U(t,*). 

We  modified  the  form  of  the  local  time  error  (per  step)  estimator  in  LSODI 
so  that  an  attempt  is  made  to  obtain 

npde  NPDE 

l  lest  (t ,.)«  <  (T0L)Z(  l  IU  (t,.)lf  +  1).  (6) 

i=l  u  i=l  u 

Here,  est(t,*)  =  {est^C  t ,  •) }  £  S  ( 5 )  is  the  function  whose  nodal 

values  are  the  local  time  error  estimates  computed  in  LSODI. 

The  character  of  the  ODE  intial  value  problem  (5)  is  assumed  to  be  such 

that 

npde  i  ,  o  npde  , 

l  9U  (t,0  -  u  (t,.)»n  <  C(T0L)^(  l  bu  (t,*)ni  +  l)  (7) 

i-1  U  i-1  U 

with  a  reasonable  constant  C  independent  of  5.  Moreover,  TOL  is  fssumed  to 
be  selected  such  that  | | |  U-U| | |  is  small  with  respect  to  | | | u-U | | | , 
and  hence  the  total  error  e  =  u  -  U  satisfies 

I  I |el I  I  =  || | u-U | | |  =  0(max  h^)  as  max  hn  +  0.  (8) 

n  n 

The  question  of  how  small  TOL  must  be  is  not  addressed  here.  For  some  results 
in  this  direction  see  Babuska  and  Luskin  [2]. 
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3.  THE  ADAPTIVE  FEMOL 

The  MOL  variation  presented  here  differs  from  the  MOL  just  described  in 
two  basic  ways.  First ,  reliable  error  information  is  obtained  during  problem 
integration  by  computing  an  a  posteriori  estimate  E(t)  of  |||e(t,*)||| 
at  each  t  in  a  set  of  initially  provided,  equally  spaced  times  {ck^k*l  K  c 
(0,T).  E  is  described  in  the  next  section.  Second,  adaptive  regridding 
occurs  at  times  {Tm)m>1  c  {t^}^  R.  The  selection  of  the  regridding  times 
depends  on  the  computed  values  of  E( • )  . 

Integration  begins  and  procedes  as  in  the  nonadaptive  FEMOL  until 
regridding  occurs  at  some  time  t^.  A  new  grid  6+  is  created  from  the 
present  grid  6  at  t^  by  uniformly  subdividing  some  elements  and  removing 
groups  of  contiguous  nodes  to  coalesce  others.  New  initial  data  U(t*,*)  ( 
S(<5+)  is  determined  simply  by  interpolating  the  already  computed 
U( t^ , • )  €  S(6)  (for  the  present  class  of  problems,  linear  interpolation  is 
more  efficient  and  has  been  observed  to  be  no  less  accurate  than  L2  or  other 
"global"  projections). 

The  grid  6  and  data  U(t^,*)  define  a  semidiscrete  approximate 
U(t,*)  €  S(6+)  for  t  >  t^  as  in  Eqs.  (4)  and  an  ODE  initial  value  problem 
of  the  form  (5)  for  solution  values  at  the  nodes  of  6+.  Integration  of  this 
problem  commences  at  time  t^  by  restarting  LSODI  with  the  same  error 
tolerance  TOL.  ODE  integration  procedes  smoothly  through  times  when 
regridding  does  not  occur.  It  continues  forward  in  time  until  regridding 
again  occurs  or  the  final  time  T  is  reached. 

Let  us  summarize  the  input  and  output  for  this  procedure.  The  input 
provided  at  the  initial  time  Tq  -  0  consists  of 
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*  a  grid  <5q, 

*  a  time  discretization  error  tolerance  TOL  >  0, 

*  a  space  discretization  error  tolerance  EPS  >  0,  and 


*  the  times  {tk}ki.1>K"  {kT/K*k«l,K  at  whlch  (*)  is  computed. 


The  output  from  the  procedure  consists  of 

*  spatial  accuracy  estimates  {E  ( t^  )  }  k_^  Rf 

*  regridding  times  (T^j  c  {tk}k=lfK> 

*  corresponding  grids  { 6m } ,  ,  and 


>  (10) 


*  an  approximate  solution  U,  which  is  in  S ( for  t  f  iTm'^m+l^'J 

There  typically  are  a  great  many  ODE  integration  time  steps  taken  on  each 

(tk,tk+j),  and  many  of  these  intervals  in  each  (Tm»Tm+l^*  ^he  T0L_ 

dependent  ODE  time  stepsize  sequence  generated  by  LSODI  generally  increases  on 

each  (^m>^m+i)»  and  It  increases  most  rapidly  just  after  being  restarted 

at  Tm»  when  LSODI  is  permitted  to  choose  the  Initial  stepsize.  As  in  the 

nonadaptive  FEMOL,  we  assume  that  the  errors  due  to  the  grids  (6  }  _  . 

m  m>  l 

dominate  those  due  to  time  discretization. 

The  primary  goal  in  selecting  (Tra}m>i  and  constructing  is  to 

obtain  a  reliable  solution  in  the  sense  that 


| I |e(t ,•)  1 1  |  <  EPS  | | | u< t , • ) M | 


t  €  (0,T) 


(ID 


The  secondary  goal  is  that  {Tm}  and  {6^}  should  be  chosen  so  that 

the  total  work  required  to  integrate  the  IBVP  on  Q  x  (0,T) 

is  less  than  that  required  with  the  same  input  and  any  other  (12) 

selections  of  regridding  times  {Tg}g>j  c  -C R  and  grids 

{6s}s>l  whlch  yield  (11)* 

The  basic  strategy  used  to  achieve  these  goals  is  summarized  as  follows: 

(i)  Integration  precedes  only  in  the  positive  time  direction  and  no 
information  is  obtained  from  the  future. 

(ii)  Information  is  collected  only  at  the  input  times  {t^}.  The  amount  of 
information  stored  at  any  time  is  small  and  is  discarded  after  one 
regridding  takes  place  in  the  future. 

(lii)  Regridding  occurs  at  some  tk  if,  and  only  if  E(tk)  >  *95 
EPS  |||U(tk,.)|||. 

(iv)  Regridding  is  carried  out  at  tk  in  order  that  E(tfc)  i  EPSDN 

I MuCtfc**) I  1 1 »  where  EPSDN  €  [.6  EPS,  .9  EPS]  is  adaptively  chosen 
and  E(tk)  is  a  prediction  (described  in  the  next  section)  of  the 
I  I) • I  I i —error  immediately  after  time  tk  with  the  new  grid. 

Pure  relative  | I  I  *  I | | —error  control  is  used  so  that  the  reader  can 
better  compare  the  method's  ability  to  estimate  and  control  errors  in  the  four 
problems  of  Section  8.  In  practice,  EPS  might  better  be  a  combined 
relative/absolute  error  tolerance,  as  the  time  discretization  error  tolerance 
TOL  is  here.  Note  that  the  constant  .95  in  (iii)  introduces  a  high  risk  of 
violating  (11)  on  some  time  interval  (tic»tk+l^*  This  constant  was  chosen  in 
order  to  show  the  dependence  of  error  control  on  error  estimation  in  the 
experiments  described  In  Section  8. 


From  (iii)  and  (iv)  one  sees  that  the  regridding  times  are  implicitly 
determined  by  EPSDN  and  the  constructed  grids.  A  grid  is  retained  until 
E/|||U|||  has  grown  to  at  least  G  times  its  value  just  after  grid 
construction,  where  G  -  .95  EPS/ EPSDN  €  (1.05,  1.60).  Let  us  remark  that  the 
character  of  the  PDEs  is  assumed  to  be  sufficiently  dissipative  so  that  the 
Influence  of  local  errors  decreases  in  time. 

In  attempting  to  achieve  the  secondary  goal  (12),  a  work  expression  is 
employed  which  utilizes  ODE  time  stepsize  information  and  computed  values  of 
£(*).  This  expression  is  used  to  select  EPSDN  and  is  described  in  Section  6. 
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4.  ERROR  ESTIMATION 

Let  t  be  one  of  the  given  input  times  (tvJk-l.K  an<*  ®  of  c^e  form 
(2),  (3)  be  the  last  space  grid  constructed  before  time  t.  The  a  posteriori 
estimate  E(t)  of  |||e(t,*)|||  *  |||u(t,»)  -  U(t,*)|||  is 

N  NPDE  .  Vo 

E(t)  -  {  I  [  |n*(t>|2}  ,  (13) 

n-1  i-1 

where  n*  is  the  local  error  indicator  for  the  ith  PDE  on  the  element: 


.x  ,+x 

0;  if  .‘(-Ei-i) 


0, 


1.2 


x 

n  i  2 

J  )r  (t,x)|  dx;  otherwise. 


(14) 


The  function  ri(t,»)  is  the  residual  of  the  ith  PDE  (neglecting  the 
discontinuities  at  the  nodes): 


rt(t,x)  -  U^(t,x)  -  a*(x)  l/(t,x)  -  f1  ( t  ,x,U(  t  ,x) )  .  (15) 

The  residual  is  obtained  via  nodal  values  in  U[t]  and  ^  U[t]  from  LSODI 
and  integrated  in  (14)  with  2-point  Gaussian  quadrature. 

The  quality  of  the  estimator  E  can  be  measured  with  the  effectivity 
index 

8(t)  =  E( t)/ | | | e(t  ,•) | 1 1 .  (16) 

This  quality  was  theoretically  analyzed  in  the  setting  of  linear  uncoupled 
PDEs  by  the  authors  [6],  [7],  where  it  was  shown  that 
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max  0(t)  1  (17) 

t 

as  the  space  grid  sizes  converge  to  zero.  The  key  assumptions  used  in  the 
proof  of  (17)  were  that  the  exact  solution  u  is  sufficiently  smooth,  u^ 
does  not  degenerate  to  the  zero  function,  ODE  integrations  are  sufficiently 
accurate,  and  that  the  grids  are  not  too  irregular  or  modified  too  frequently. 

All  theoretical  details  have  not  been  carried  out  for  nonlinear  reaction- 
diffusion  systems,  but  computational  experience  suggests  that  similar  theory 
applies  here.  The  evaluation  of  E  is  included  in  many  experiments  in 
Section  8. 

Let  us  now  assume  that  t  €  {Tm}_  .  and  the  grid  6  is  about  to  be 

m  m>  1 

modified.  The  | | | • | | | -error  immediately  following  a  proposed  regridding 


6+  -  On"  -  x*  <  x|  <  •••  <  x+  -  3fl+} 

1  N 


is  estimated  with 


E(t+) 


N+  NPDE 

{  l  l 

j-1  i-1 


(t+)|2fy2 


» 


(18) 


where  {hj(t+))  are  predicted  values  of  the  error  indicators  for  the  grid 
6+«  These  indicators  are  determined  from  the  already  computed  (q*(t)}  in 
the  following  way. 

If  some  element  (xn_i»xn^  is  t0  ^  refined  into  the  union 


U 

J“J0 


+ 


) 


of  q  uniform  subelements,  then  a  q-fold  decrease  in  the  contribution  to  the 
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[•Ill-error  on  (xn_i»xn)  is  predicted: 


V1+q  NPDE 


I  l  !nj(t+)|2 


NPDE 


j-jQ  1-1 


q'2  l  h*(t) | 2, 

i-1  n 


(19) 


Alternatively,  if  some  q  contiguous  elements  {(x„  , ,x„)} 

1  n— l *  n  n— np, 

present  grid  are  to  be  coalesced  to  form  one  new  element  (x^  j. 


„  in  the 

no  i+q 

Xj),  then 


NPDE 

l 


i-1 


Inj(t+)|2 


no‘1+q 

(  I 


no~1+<* 

l 

n-n^ 


NPDE 


l 

i-1 


\\(t) 


2 


(20) 


If  it  happens  that  h^  -  *••  -  then  (20)  corresponds  to  a 

predicted  q-fold  increase  in  the  contribution  to  the  | | | • | | (-error  on 
(xnQ-l»xn0-l+q)  *  (xj-l»x+p* 

These  predictions,  as  those  used  in  [7],  are  based  on  the  expectation  that 
for  n  -  1  ,N 


NPDE 

l  lnn(t) 

i-1 


2 


h2  Xn  NPDE  1 

*  12  /  l  a  ^x^!uxx^t,x)'  dx  *  0+°(hn))  as  hn  ♦  0.  (21) 

Vi  1-1 

When  the  grid  6  is  to  be  modified,  additional  information  is  extracted 
from  the  error  indicators  (n*(t)}.  a  piecewise  constant  function  is 
constructed : 


NPDE 


1  *  wU  •  -  .  i. 

5  h» ,L  iv*>n  •  ■  ( “-i.*. 


w(t ,x) 


(22) 


which  Is  an  approximation  of  the  cube  root  of  the  integrand  in  (21).  From  the 
definition  (13)  of  E( t),  one  sees  that 

„  l2  x 

.  N  h  n 

EZ(t)  -  l  ~  /  w3(t,x)dx.  (23) 

n’1  xn-l 

The  function  w(t,»)  and  (23)  are  explicitly  used  in  determining  the  "shape" 
of  the  new  grid  6+ • 
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5.  GRID  SHAPE  AND  INTENSITY 

In  refining  and  derefining  a  grid  6,  two  properties  of  5  are 
changed:  its  shape  and  its  intensity.  To  explain  what  is  meant,  we  begin 
more  generally  by  defining  the  shape  of  a  positive  integrable  function  £  on 
ft  as  the  graph  of  the  function 

R(x,y)  =  £(x)/£( y) ;  x.y  €  ft 

and  the  intensity  of  £  simply  as 

I  =  /  £(x)  dx. 

ft 


Two  positive  functions  have  identical  shape  if  and  only  if  they  are  constant 
multiples  of  one  another.  A  positive  integrable  £  can  be  magnified 
(multiplied  by  a  constant)  to  get  a  function  having  any  positive  intensity  and 
the  shape  of  £. 

The  shape  and  intensity  of  a  grid  5  are  taken  to  be  those  of  the 
associated  grid  function 


VX) 


'l/h  :  x  6  (x  ,x  )  for  n  »  1,N, 
n  n—  l  n 

„.5(l/h  +  l/h  );  x  -  x  for  n  -  1,  N  -  1. 

n  n+l  n 


(24) 


It  is  easily  checked  that  6  has  intensity  E^  »  N. 

Let  us  briefly  examine  the  effects  of  grid  shape  and  intensity  on  the 
estimated  | | | • | | | —error  in  the  present  method.  To  this  end,  at  a  time  t 
when  a  grid  5  is  being  used,  we  may  use  £(g  in  the  expression  (23)  for 
E^(t)  to  obtain 

E(t)  -  A(w(t,.),56) 
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w3(t ,x) 
C^(x) 


(25) 


Since  w(t,»)  (cf.  (22))  tends  toward  a  grid-independent  function  as  the 
information  defining  it  improves,  let  us  consider  it  to  be  a  known  function. 

Now,  had  each  element  of  5  been  bisected  at  the  grid's  creation,  for 
example ,  the  resulting  grid  would  have  had  the  same  shape  and  twice  the 
intensity  of  6.  One  sees  by  examining  the  integrand  in  (25)  that  E(t) 
would  have  been  half  as  large.  Alternatively,  if  the  N  -  1  interior  nodes 
of  6  had  been  rearranged — i.e.  the  same  number  of  extra  nodes  added  as 
deleted  here — then  the  shape  of  6  would  have  been  different,  but  the 
intensity  would  not.  The  change  in  E(t)  then  would  have  depended  strongly 
on  how  the  shape  of  6  had  changed  with  respect  to  the  shape  of  w( t , • ) . 

In  constructing  a  grid,  all  reasonable  refinement/derefinement  algorithms 
for  time-dependent  PDEs  must  select  a  shape  and  an  intensity  for  the  grid 
which  will  work  in  the  future.  This  is  usually  done  indirectly — by  focussing 
completely  on  local  error  estimates,  for  example.  In  the  adaptive  FEMOL, 
these  two  properties  are  controlled  more  directly.  A  grid  is  constructed  by 

*  explicitly  selecting  a  positive  model  grid  function  £, 

*  magnifying  £  to  yield  an  implicitly  defined  model  grid  intensity  I, 
and 

*  refining  and  derefining  so  that  the  resulting  grid  has  a  shape  closely 
resembling  that  of  £,  with  the  number  of  elements  approximately  equal 
to  I . 

Let  Tm  be  the  time  at  which  a  grid  5+  is  to  be  constructed  from  6. 
Recalling  the  basic  regridding  strategy  described  in  Section  3,  we  know  that 


■  '♦ 
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E(  Ta) 

ll|U(Tm.*>|||  > 


.95  EPS 


and  that  the  predicted  error  immediately  following  time  Tm  with  6+  is  to 
satisfy 


in 


l!u(Tm.-)||| 


EPSDN, 


(26) 


with  EPSDN  €  [.6  EPS,  .9  EPS]  still  to  be  chosen.  The  model  intensity  I  is 
defined  by  (26).  I  is  the  minimum  intensity  that  6+  can  have  while  having 
the  shape  of  £  and  yielding  (26). 

Assume  that  £  and  EPSDN  have  already  been  selected.  The  construction  of 
6+  is  then  carried  out  as  follows.  The  elements  ( (xn-l  ,xn^n=l  N  of  5 


are  uniformly  subdivided  or  coalesced  according  to  desired  new  local  element 

f 

2 


x  +x 

sizes  at  the  midpoints  {xn_  y  }  =  {— n~  — -}  of  6: 


h+(x 


n-  V?  ^ 


1 


C^(xn-lA>  ’ 


n  =>  1 ,  N. 


(27) 


The  constant  C  »  C(EPSDN)  >0  in  (27)  is  initialized  in  a  predetermined  way 
and  iteratively  adjusted  until  (26)  holds.  At  each  stage  of  the 
iteration,  E(T*)  is  obtained  from  the  already  computed  local  indicators 
(.^(T^)}  for  6  as  described  in  Section  4,  with  refinement  indices  gotten 
from  the  ratios  {hn/h+(xn_  \j  )}  of  old  to  new  element  sizes. 

Further  details  of  this  iteration  algorithm  are  not  given  here.  It  is 
noted,  however,  that  none  of  the  usual  mesh  constraints  (on  element 
subdivision  ratio,  ratio  of  neighboring  element  sizes,  number  of  contiguous 
elements  coalesced,  etc.)  are  explicitly  imposed.  Some  are  incorporated 
naturally  through  the  selection  of  £. 
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In  view  of  (27),  one  can  expect  that  the  grid  function  £  (•)  associated 

+  « 
with  6  will  resemble  the  function  <>£(•)  in  shape  and  intensity. 

Furthermore,  one  expects  as  in  (25)  that 


E(t)  L  A(w(t,.),C£) 

*  77  A(w(t, •),£);  t  >  T  , 

u  m 

and  from  (26)  that 


(28) 


1 

c  tniu(T.,.)iir 


EPSON. 


(29) 


From  (28)  it  is  clear  that  £  should  reflect  the  present  and  predicted 
future  shape  of  w(t,*).  £  is  chosen  to  be  a  majorant  of  w(Tm,»)  in  the 

adaptive  FEMOL.  The  selection  process  anticipates  the  future  through 
comparison  of  a  small  amount  of  information  collected  at  and  the  previous 

regridding  time  Tm_| .  This  process  is  described  in  Section  7. 

The  intensity,  or  number  of  elements  in  5+  (and  number  of  ODEs  in  time 
to  be  solved)  depends  both  on  £  and  the  final  value  of  C  in  (27).  Once  £ 
has  been  selected,  however,  the  intensity  is  determined  by  the  parameter 
EPSON,  as  evidenced  by  the  inverse  relation  between  C  and  EPSDN  in  (29). 

The  selection  of  EPSDN  is  now  described. 
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6.  CONTROL  OF  GRID  INTENSITY 

Assume  that  the  model  grid  function  5  has  been  chosen  at  a  time  Tffl 

when  a  grid  <S+  is  to  be  constructed  from  6.  The  parameter  EPSDN 

determining  the  intensity  of  6  is  selected  In  an  attempt  to  minimize  the 

work  per  future  (unknown)  time  step  Tm+^  -  which  will  be  required  with  a 

grid  having  the  shape  of 

Consider  two  extreme  values  of  EPSDN.  If  EPSDN  is  chosen  to  be  .6  EPS 

(smallest  permitted  value),  ^m+1  ~  ^m  would  be  as  large  as  possible,  but  so 

also  would  be  the  number  of  ODEs ,  since  this  number  depends  inversely  on 

EPSDN.  With  EPSDN  *  .9  EPS  (largest  permitted  value),  the  smallest  ODE  system 

would  be  integrated  for  the  shortest  period  of  time.  The  sum  of  such 

integrations  could  be  quite  costly,  not  only  because  of  grid  construction 

costs,  but  also  since  efficient  use  of  an  ODE  solver  whose  internal  stepsizes 

increase  is  made  only  if  the  solver  is  not  frequently  restarted. 

If  Tm  ■  Tj  (first  regridding) ,  EPSDN  is  taken  to  be  .9  EPS.  Otherwise, 

it  is  selected  by  estimating  the  marginal  benefit  which  would  have  resulted, 

had  the  value  EPSDN_  ,  of  EPSDN  used  at  T  ,  been  altered.  Let 
tn—  l  m—  l 

E(t)  -  ( - — 5  ---S - )  •  EPSDN 

sup  E(s)/| | |U( s , • ) | | |  m 

s€(Tm-l»t) 

for  t  >  Tm_j .  The  number  c(t)  is  an  estimate  of  the  largest  value  which 

EPSDN  could  have  taken  at  while  leading  to  successful  PDE  integration 

(i.e.  E/ 1  |  |  U 1 1 1  <  .95  EPS)  on  (Tm_pt).  Note  that  e(*)  decreases  on 

(Tm_i,Tm)  and,  by  the  basic  regridding  strategy  described  in  Section  3,  that 

e(T+  ,)  1  .95  EPS  and  e(T  )  <  EPSDN- 
v  ra-i '  in  m- 1 

The  work  per  regridding  time  step  which  would  have  been  required  by 
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selecting  EPSDN  -  e(t)  at  Tm_j  can  be  expressed  as 


WORK(t) 


(STEP(t)  +  c) 


Here,  STEP(t)  is  the  number  of  ODE  time  steps  taken  on  (Tm_j,t)  with  EPSDN  » 

EPSDN  The  constant  c  is  related  to  the  overhead  incurred  at  T  .  in 

m-  i  n»- 1 

grid  construction  and  data  initialization. 

Using  stepsize  information  returned  from  the  ODE  solver,  computed  values 

of  £(•)  and  simple  extrapolation  for  t  >  Tm ,  profiles  of  £(•)  and 

WORK(*)  are  updated  as  time  increases  and  are  used  to  determine  an  "optimal" 

value  of  EPSDN  at  T  ,  in  retrospect.  This  value  then  governs  the  selection 

m- 1 

of  EPSDN  at  T  .  Details  of  the  final  selection  are  not  important.  What  is 

m 

important  is  that  changes  in  EPSDN  smoothly  reflect  observed  trends.  Here, 
EPSDN  is  increased  (decreased)  by  no  more  than  .02  EPS  (.05  EPS)  from  one 
regridding  time  to  the  next. 
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7.  CONTROL  OF  GRID  SHAPE 

As  indicated  in  Section  5,  the  shape  of  the  grid  6+  constructed  from  6 
at  time  Tm  resembles  that  of  the  model  grid  function  £,  which  is  chosen  as 
a  majorant  of  the  function  w(Tm,*)«  Recall  that  w(Tm,*)  is  an 
approximation  of  a  function  related  to  the  exact  solution's  second  space 
derivatives  at  Tm  (cf.  (22))  whose  values  are  extracted  from  the  local  error 
indicators  for  6.  The  process  of  selecting  ?  is  introduced  by  considering 
two  extreme  possibilities: 

(i)  ?(•)  -  v(Tm.O 

and 

(ii)  £(•)  =  max  w(T  ,x)  (i.e.  a  constant  function). 

x*n  m 

Consider  £  as  in  (i) .  Using  the  representation  (28)  for  the  predicted 
error  estimate  E(T*)  immediately  after  Tm  and  standard  arguments  of  the 
calculus  of  variations,  as  in  Babu§ka  and  Rheinboldt  [3],  it  can  be  shown  that 
the  resulting  grid  <5+  would  have  the  least  intensity  of  all  grids  lowering 
the  estimated  relative  ) | | *  j | | —error  to  any  given  EPSDN  value.  Since  nothing 
is  known  for  certain  about  the  future  shape  of  w(t,»)  and  this  choice  yields 
the  ODE  system  of  smallest  possible  size,  one  might  call  (i)  the  low  cost 
alternative . 

The  local  error  indicators  would  be  approximately  equilibrated  with  (i), 
but  only  near  time  Tm  if  the  shape  of  w(t,»)  rapidly  changes.  Because  an 
immediate  and  costly  regridding  could  be  necessitated,  the  choice  (i)  is  not 
"optimal"  in  any  practically  useful  sense  of  the  word  for  parabolic  PDEs  whose 
solutions'  higher  derivatives  undergo  changes  in  shape.  For  such  problems, 

(i)  might  be  better  termed  the  high  risk  alternative.  The  local  orientation 
of  (i)  introduces  a  strong  dependence  of  error  control  on  the  input  parameters 
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(i.e.g.  the  probability  of  a  failure  occuring  on  some  interval  i3 

high)  and  a  strong  dependence  on  each  piece  of  local  information  (it  is  often 
Che  case  that  a  "few"  of  the  local  indicators  defining  w(Tm»*)  are 
relatively  very  inaccurate) . 

The  choice  (ii)  of  £  represents  the  opposite  end  of  the  risk-cost 
spectrum.  It  is  a  high  cost  alternative,  since  the  resulting  grid  6+  would 
be  nearly  uniform  and  have  the  largest  intensity  of  all  reasonable  grids 
lowering  E/|))u||)  to  any  given  EPSDN  value.  It  imposes  a  low  risk  of 
losing  control  over  errors,  since  unpredictable  changes  in  the  shape  of  higher 
solution  derivatives  would  be  accounted  for,  as  well  as  isolated  instances  of 
inaccurate  local  information. 

Neither  (i)  nor  (ii)  represents  a  viable  alternative.  Heuristics  must  be 
employed  to  balance  risk  and  cost  and  to  predict  the  future.  All  reasonable 
algorithms  of  the  present  type  must  use  heuristics  to  predict  the  future, 
since  the  future  is  generally  far  ahead  of  the  present — much  farther  than  for 
which  accurate  and  inexpensive  local  extrapolations  can  be  used. 

Physical  reasoning  can  be  used — i.e.  mass  conservation,  expected  wave 
speeds,  etc. — or  numerical  simulations  can  be  adopted — i.e.  global  extrapola¬ 
tions,  averaged  or  extended  information,  etc.  Can  such  prediction  processes 
ever  be  quantitatively  assessed  or  compared,  or  better  justified  in  a  general 
way?  If  so,  the  first  step  may  be  to  see  them  for  what  they  are — pattern 
recognition  processes.  Such  a  process  consists  of  three  (generally 
nondistinct)  stages: 

*  Representation  -  reduction  of  (perhaps  "noisy")  data  into  a  convenient 
and  invariant  form, 

*  Feature  Extraction  -  relevant  measurements  taken  from  the  reduced 


data,  and 
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Classification  -  decisions  made  by  comparing  feature  values  in  an 
attempt  to  improve  recognition  or  to  avoid  misrecognition . 


Let  us  use  this  framework  and  terminology  to  describe  the  construction  of 
£  at  time  T  .  The  "pattern"  we  wish  to  predict,  or  recognize  is  the  shape 
of  the  function 


w(x) 


max 
te [T  ,T 


m  nH-1 


1 


w( t,x)  , 


where  Tm+j  is  unknown ,  but  Tm+^  -  Tm  is  assumed  to  be  comparable  to  Tffl  - 
Tm_p  The  attempt  to  do  this  consists  of  constructing  a  piecewise  linear  £ 
which  majorizes  w(Tm,*)  and  whose  shape  approximates  that  of  w.  The 
reasoning  behind  this  choice  comes  from  a  variational  principle  related  to  the 
present  problem  and  which  involves  the  functional  A  used  in  the 
representation  (28)  of  E. 

In  order  to  construct  £,  each  finite  element  grid  Is  required  to  contain 
the  fixed  nodes  of  a  uniform  macro  grid  A.  A  is  supplied  as  input  at  the 
initial  time  and  is  a  3- level  grid,  whose  several  "large"  level  1  macro 
elements  each  have  size  4H,  and  whose  level  2  and  3  subelements  have  sizes  2H 
and  H,  respectively.  The  size  H  (4H)  is  related  to  the  maximum  (minimum) 
risk  of  losing  control  over  ||(e||j  that  one  is  a  priori  willing  to  take  and 
the  minimum  (maximum)  price  one  Is  willing  to  pay  to  keep  it.  The  algorithm 
for  constructing  £  tries  to  manage  risk  and  cost  on  each  level  1  macro 
element  (X,X+4H)  by  first 

reducing  data  -  The  many  pieces  of  data  defining  w(Tm,»)  on 
(X,X+4H)  are  replaced  by  three  piecewise  constant  functions 

2  3*  w^ere  W^(Tm>»)  takes  the  maximum  value  of  w(Tm,»)  on 
each  of  the  2U  *  level  u  macro  elements  contained  in  (X,X+4H).  In  a 


I 


similar  manner,  three  piecewise  constant  representations 

(W  (T  ,•)}  ...  of  w(T  .,•)  on  (X,X+4H)  were  formed  at  time 

y  m-1  y«l,2,3  m-1’ 

T  ,,  where  w(T  .,•)  =0  if  T  ,  *  0  (the  initial  time). 

m-i  ra-i  ui“i 

What  relevant  features  can  be  extracted  from  these  representations?  In 
solving  many  parabolic  reaction-diffusion  problems,  it  has  been  observed  that 
there  often  is  a  correlation,  on  some  scale,  between  spatial  differences  in 
w(Tm_p«)  and  the  way  w  subsequently  grows  in  time.  The  algorithm  looks 
for  such  a  correlation  on  (X.X+4H)  x  (Tm-l»^m^  by  taking  three  measurements: 


X+4H 

/  i«3<v*>  -  vvi-x>idx 


y  *  1.2,3. 


These  three  feature  values  are  used  to  classify  the  evolution  of  w  on 
(X,X+4H)  x  (Tm_i»Tm).  An  "active"  level  is  taken  to  be  that  corresponding  to 
the  largest  index  y*  (  {1,2,3}  for  which 


M  *  <  M  ; 

y*  y 


y  -  1,2,3. 


If  y*  *  3  (which  always  is  the  case  at  the  first  regridding  time  T^), 
either  w  did  not  grow  on  (X,X+4H)  x  (Tffl_^ ,Tm)  or  it  is  concluded  that 
spatial  differences  in  w  at  Tm_j  were  not  the  source  of  its  growth.  If 
y*  *  2  (1),  it  is  concluded  that  either  a  clear  correlation  existed  on  the 
scale  2H  (4H)  or  w  evolved  in  a  way  which  was  unpredictable  at  Tm_^  on 
any  smaller  space  scale. 

Determining  an  "active"  level  y*  for  (X.X+4H)  can  be  interpreted  as 
choosing  one  of  three  directions  in  the  (x,t)-plane.  The  directions  are 
defined  by  the  magnitudes  of  time  and  space  differences, 

(W3(Tra’*)  ~  W3^Tm-l’*^  and  ^W3^Tm-l’*^  ”  Wy^Tra-l’*^  of  W3’  which  is  the 
most  local  of  the  three  representations  of  w.  The  direction  corresponding  to 
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U*  is  that  along  which  Wg(t,*)  changed  the  least  for  t  6 

Having  classified  the  past  evolution  of  w  on  each  level  1  macro  element 
in  A ,  these  classifications  are  used  to  predict  the  extent  to  which  spatial 
differences  in  w(Tm,*)  will  affect  the  size  of  w(t,»)  in  neighboring 
regions  for  t  >  "1^.  A  macro  subgrid  consisting  of  the  boundary  nodes  of 
every  large  level  1  element  and  the  boundary  nodes  of  every  "active"  level  2 
or  3  element  is  formed  (i.e.g.  if  \i*  *  2  for  some  (X.X+4H),  then  the 
node  X  +  2H  is  in  the  subgrid) .  The  model  grid  function  £  is  taken  to  be 
the  piecewise  linear  function  on  the  subgrid  whose  jc^  subgrid  nodal  value  is 
equal  to  the  maximum  of  between  the  j-lst  and  j+l3t  subgrid  nodes. 

These  values  ari  obtained  from  those  of  the  piecewise  constant  representations 
{W  (Tffl,»)}>  which  are  subsequently  stored  for  use  in  computing  feature 
values  at  the  next  regridding  time. 

There  are  two  important  effects  of  choosing  5  to  be  continuous  and 
piecewise  linear,  as  opposed  to  being  piecewise  constant.  First,  the 
dependence  of  the  selection  process  on  the  input  macro  grid  A  is  lessened. 
Second,  the  finite  element  grid  6  is  sraoothened,  since  the  local  element 

size  h+(x)  at  x  €  SI  is  inversely  proportional  to  £(x)  . 

We  show  how  the  construction  of  the  model  grid  function  £  is  carried  out 
in  two  ways.  Figures  l,  2  and  3  picture  the  results  of  simulations,  in 
which  w  was  an  exactly  known  function  related  to  the  unimodal  Gaussian 
probability  density  function.  Figures  4  and  5  picture  the  construction  of  £ 
and  the  resulting  finite  element  grid  6+  in  two  of  the  examples  in  Section 
8.  The  function  w  in  these  examples  was  extracted  from  the  local  error 
indicators,  as  described  in  Section  4.  In  each  of  Figures  1-5,  A  is  the  3- 
level  macro  grid  whose  level  1  elements  are  delineated  by  vertical  hash  marks. 
In  the  simulations  pictured  in  Figures  la  and  lb,  w  shifts  to  the  right 
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by  H  and  2H  units,  respectively,  much  as  the  second  derivative  of  a 

traveling  wave  solution  might.  In  both  cases,  £  localizes  the  left  side 

of  w(Tm,»)  as  much  as  possible,  since  w  has  not  grown  in  that  region.  The 

algorithm  "sees"  correlations  (on  different  scales)  between  spatial 

differences  in  w(T  .,•)  and  the  growth  of  w  at  the  right,  however,  and 

m-1 

the  extensions  of  £  bo  the  right  show  how  the  movement  of  w  is  predicted. 

In  a  second  pair  of  simulations ,  w  was  taken  to  be  a  skewed  Gaussian 
function  whose  width  grows  from  time  Tn_j  to  the  present  time  T  .  Figures 
2a  and  2b  show  how  the  model  grid  function  £  was  constructed  when  the  width 
of  w  doubled  and  tripled,  respectively.  As  in  the  previous  simulations,  the 
scale  of  the  correlation  "seen"  by  the  algorithm  and  used  to  predict  the 
future  shape  of  w  depends  on  how  much  w  has  grown  in  the  past . 

In  the  simulations  pictured  in  Figures  3a  and  3b,  the  intensity  of  w 
grew  by  25%  and  100%,  respectively,  but  the  shape  of  w  did  not  change.  The 
shape  of  £  closely  resembles  that  of  w  in  the  first  simulation,  but  much 
less  so  in  the  latter.  There,  the  influence  of  the  peak  of  w(Tmi*)  is 
predicted  to  spread  a  distance  4H  before  the  next  regridding  time,  since  the 
information  available  to  the  algorithm  at  Tm_j  was  inadequate  to  predict  the 
growth  of  w  on  any  smaller  space  scale .  One  could  tune  the  algorithm  so 
that  £  would  have  a  near  "optimal"  shape  for  this  type  of  problem,  but  not 
without  lowering  the  algorithm's  performance  in  more  interesting  and  complex 
problems. 

We  now  turn  to  two  examples  taken  from  Section  8.  Figure  4  shows  the 
model  grid  function  £  and  new  finite  element  grid  6+  for  Example  4,  in 
which  temperature  and  a  single  species  concentration  propagate  in  a  flame 
front  from  right  to  left.  Here,  <5+  is  the  grid  pictured  at  Tm  ■  3.8  * 

10”^  in  Figure  14  of  Section  8. 
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By  looking  ahead  of  the  moving  front,  where  w  is  predicted  to  grow  in 
the  future,  we  see  how  the  shape  of  £  determined  the  shape  of  the  grid  5*. 
While  £  is  constant  on  the  second  level  1  macro  element  pictured,  the 
corresponding  nodes  of  6+  are  not  quite  uniformly  spaced,  since  regridding 
consists  of  refining  and  derefining  the  (not  pictured)  previous  grid.  The 
parameter  EPSDN  determining  the  intensity  of  <5+  was  chosen  to  be  .85  EPS  in 
the  situation  pictured  in  Figure  4.  Had  the  relative  | | | • | | | —error  been 
reduced  to  approximately  .6  EPS  by  taking  EPSDN  equal  to  this  value,  for 
example,  3/2  as  many  nodes  would  have  been  distributed  in  the  same  relative 
way  as  pictured.  The  majority  of  the  extra  nodes  would  have  been  located  _in^ 
the  front  at_  Tm,  where  they  never  will  be  needed  to  lower  the 
111*111 -error .  The  remaining  portion  placed  ahead  of  the  front  would  only 
extend  the  next  regridding  time  by  a  small  amount. 

Figure  5  shows  the  function  £  and  the  grid  6+  which  were  constructed 
in  the  population  ecology  model  In  Example  1  of  Section  8.  6+  is  the  grid 
pictured  in  Figure  7  at  T  ■  .6,  when  the  populations  are  in  the  midst  of 
their  evolution  from  one  localized  spike  to  a  spatially  oscillatory  steady 
state  with  seven  relative  maxima  spread  over  the  entire  domain. 

The  "active"  levels  in  the  3-level  macro  grid  A  were  set  equal  to  2  in 

the  two  central  level  1  elements  at  time  T„  and  1  in  all  others.  This 

m 

reflects  the  fact  that  the  spatial  scale  on  which  the  growth  of  w  is 
correlated  with  differences  of  w(Tm_j,»)  is  larger  on  the  sides  of  the  domain 
than  at  the  center.  As  in  Figure  4,  we  see  how  the  shape  of  £  determined 
the  shape  of  6+.  The  parameter  EPSDN  determining  the  intensity  of  was 

chosen  to  be  .85  EPS  at  Tm,  as  it  was  in  the  previous  example.  Unlike  in 
the  propagating  flame  problem,  however,  where  EPSDN  was  changed  very  little 
from  one  regridding  time  to  the  next,  EPSDN  was  steadily  decreased  in 
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regrlddings  subsequent  to  the  one  pictured  in  Figure  5.  It  becomes 
computationally  more  efficient  to  increase  the  grids'  intensities  as  their 
shapes  become  more  constant  and  predict  the  spreading  influence  of  w. 

The  strengths  of  coupling  multi-grid  feature  extraction  with  implicit,  low 
order  prediction  in  time  lie  in  the  generality  of  the  approach,  its  empirical 
success  in  predicting  shape,  and  the  potential  application  of  similar 
techniques  in  two  or  three  space  dimensions.  It  is  not  important  that 
binary,  3-level  macro  grids  are  used.  Experience  indicates  that  the  process 
of  keeping  risk  and  cost  low  through  implicit  shape  prediction  improves  with 
p- level  macro  grids  as  p  is  increased,  provided  a  reasonable  number  of  large 
level  1  elements  are  present.  It  does  seem  important,  however,  that  explicit 
shape  predictors  using  representations  of  w  at  past  regridding  times  not  be 
used.  The  mesh  sizes  H  and  {Tm_Tm-l^m>l  are  not  asymptotically  small. 
Explicit  Taylor  expansions  using  these  mesh  sizes  can  be  very  inaccurate  and 
readily  destabilize  the  process  they  are  meant  to  control. 
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8.  COMPUTATIONAL  EXPERIMENTS 

The  results  of  applying  the  adaptive  FEMOL  to  numerically  solve  four 
reaction-diffusion  problems  are  described  here.  These  problems  are  of  varying 
difficulty  and  were  chosen  in  order  to  give  a  representative  evalution  of  the 
method's  performance. 

Most  of  the  experiments  fall  into  one  of  two  categories.  In  the  first,  a 
reasonable  set  of  input  parameters  was  fixed  in  each  problem  and  the  adaptive 
method  was  applied  with  a  decreasing  sequence  of  space  and  time  discretization 
error  tolerances  EPS  and  TOL,  where  each  TOL  value  was  small  with  respect  to 
the  corresponding  EPS  value.  In  the  experiments  in  the  second  category,  these 
same  TOL  values  were  used  with  the  same  method,  but  in  which  no  feedback 
related  to  regridding  was  processed,  a  posteriori  error  estimates  were  only 
computed  at  few  target  times,  and  increasing  numbers  of  uniform,  time- 
independent  finite  elements  were  used. 

The  most  important  component  of  the  adaptive  method  is  the  a  posteriori 
error  estimator  E(t)  (cf.  Section  4),  on  which  most  automatically  made 
decisions  are  based.  The  quality  of  E(t)  is  examined  for  both  changing  and 
unchanging  grids  with  the  effectivity  index  0(t)  »  E( t)/ j | | e( t) j | j  ,  where  a 
very  accurate  approximation  was  used  to  compute  the  error  e  *  u  -  U  when  the 
exact  solution  u  was  not  available.  Let  us  recall  from  Section  2  that 
|||*||j  is  a  weighted  L2  gradient  norm  (weighted  norm),  and  that 

|||e|||  is  predicted  by  theory  to  converge  linearly  to  zero  with  decreasing 
local  space  grid  sizes. 

The  primary  goal  (cf.  (11),  Section  3)  of  the  adaptive  FEMOL  is  to  keep 
the  relative  I  I  I • I 1 1 —error  below  the  tolerance  EPS  at  all  times.  The  method 
is  relatively  successful  in  achieving  this  objective,  and  we  shall  see  that 
achievement  of  this  goal  depends  strongly  on  the  quality  of  E.  Rough 
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estimates  of  the  cost  in  controlling  the  error  are  gotten  by  comparing  the  CPU 
times  required  in  the  adaptive  and  nonadaptive  runs.  These  times  include  a 
small  amount  required  in  I/O  associated  with  the  experiments,  but  within  each 
experiment  they  provide  a  fair  estimate  of  the  total  cost. 

The  secondary  goal  (cf.  (12),  Section  3)  of  the  adaptive  FEMOL  is  to  keep 
costs  as  low  as  possible,  given  a  set  of  input  and  the  control  mechanisms  of 
the  method.  Note  that  the  CPU  time  comparisons  mentioned  above  do  not  provide 
a  fair  measure  of  performance  for  this  goal.  At  the  end  of  this  section  we 
summarize  results  of  an  experiment  designed  to  more  adequately  gage  the 
relative  cost  of  the  adaptive  procedure. 

All  experimental  runs  were  made  using  the  program  FEM0L1,  double  precision 
arithmetic,  and  the  FORTRAN  H  compiler  on  the  IBM  System  370/3081K  at  the 
NIH.  The  notation  used  to  describe  the  experimental  results  has  been 
introduced  earlier  or  is  self-evident,  with  the  exceptions  of 

CPU  *  IBM  308 IK  seconds, 

mean  NO.  ELTS.  *  (1/T)  •  Y  N  , (T  -T 

u  ra— i  m  ra— I  ’ 

1 

where  Nm_ ^  is  the  number  of  finite  elements  used  from  time  Tffl_j  to  time 
Tra,  and  A(J,2J,4J),  which  denotes  a  3-level  macro  grid  containing  J,  2J 
and  4J  level  1,  2  and  3  macro  elements,  with  a  total  of  4J  +  1  fixed, 
uniformly  spaced  nodes  (cf.  Section  7,-. 

Example  1.  Population  Ecology  Model 

The  system  considered  here  was  proposed  by  Murray  [19]  to  model  certain 
planktonic  predator-prey  situations  in  which  crowding  is  a  factor.  The 
predators  v  (zooplankton)  and  essentially  static  prey  u  (phytoplankton) 
satisfy 
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u  -  .0125  u  ■  [f(u)  -  v]u 

t  XX 


vt  -  vxx  =  [u  -  g(v)]v 


t  >  0 

X  €  (0,2.5) 


(30) 


u  =  v  =0 

X  X 


t  >  0 

x  -  0,  2.5 


(31) 


and 


where 


u  =  u0,  V  -  v0 


f(u)  =  (35  +  16u  -  0/9 


t  *  0  (32) 

x  6  (0,2.5) 


(33) 


g(v)  =  (5  +  2v)/5 


and  the  initial  populations  {uq»vo^  are  as  Picture<*  *n  Figure  6,  along  with 
the  steady  state  to  which  {u,v}  evolves.  There  are  many  steady  state 
solutions  of  Eqs .  (30),  (31)  (the  stable  solution  {u,v}  =  {5,10}  of  the 
diffusionless  ODE  system,  for  example),  and  the  one  to  which  {u,v}  evolves 
depends  on  which  eigenfunction  of  the  linearized  coupled  elliptic  operator  the 
perturbation  of  the  initial  data  from  {5,10}  most  resembles.  These 
piecewise  linear  initial  data  were  represented  exactly  in  all  of  the 
experiments  described  here. 

This  problem  is  relatively  easy,  and  a  properly  selected  number  of 
uniformly  spaced  grid  points  would  adequately  keep  the  errors  below  a  desired 
tolerance.  The  role  of  the  adaptive  FEMOL  in  this  problem  is  that  of  a 
convenient  and  reliable  tool,  with  which  the  probability  of  success  in 
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achieving  the  goal  (11)  in  one  application  is  high. 

Tables  I  and  II  summarize  the  results  for  this  problem.  In  computing  the 
relative  | ! i *  I | 1 —error  and  9  at  the  10  target  times  for  the  entries  of 
these  tables,  the  "exact"  solution  u  =  {u,v}  was  taken  to  be  the  numerical 
solution  computed  with  512  uniform  elements  and  very  accurate  ODE  time 
integration  in  the  nonadaptive  FEMOL.  The  approximations  computed  with  the 
input  parameters  in  rows  3  and  4  of  these  tables  were  too  accurate  to  be 
compared  with  u.  The  corresponding  listed  relative  | | | “ | | | —errors  were 
therefore  taken  to  be  those  estimated  in  the  method. 

The  value  t  in  Table  I,  row  1  was  chosen  as  the  largest  value  of  TOL  for 
which  the  maximum  relative  | | | “ | | | —error  due  to  time  discretization  was  less 
than  .01.  All  other  input  TOL  values  in  Tables  I  and  II  were  chosen  by 
subdividing  t,  as  might  be  done  in  practice.  This  same  procedure  was 
repeated  in  Examples  2  and  3. 

Note  the  very  high  quality  of  the  error  estimator  E  in  this  example,  as 
seen  in  the  | 0—1 |  columns  of  both  Tales  I  and  II.  In  no  Instance  did  the 
estimated  ! | | ” l tl —error  differ  from  the  true  111*111 —error  by  over  3.4%. 
This  difference  evidently  diminishes  as  the  relative  | | | • | | |-error  +■  0,  more 
rapidly  in  the  uniform  nonadaptive  case  than  in  the  adaptive. 

We  next  observe  that  the  adaptive  FEMOL  successfully  achieved  its  primary 
goal  in  all  4  runs,  as  witnessed  by  comparing  the  desired  maximum  relative 
I  I  I  * ! II -errors  in  the  EPS  column  of  Table  II  with  the  actual  maximum  relative 
111*111 -errors •  In  doing  this,  the  mean  number  of  elements  used  in  each  run 
exceeded  the  corresponding  number  of  uniform  elements  by  5  -  25%.  The  actual 
surcharge  was  about  100%,  however,  as  is  seen  by  comparing  CPU  times  in 
Tables  I  and  II.  We  mention  that  the  selection  of  the  uniform  finite  element 
grid  providing  a  desired  accuracy  EPS  could  be  expensive,  and  this  cost  is  not 
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reflected  in  the  CPU  entries  of  Table  I.  The  adaptive  method  must  have  higher 
overhead  costs  than  the  nonadaptive  method,  since  the  adaptive  method  makes 
decisions  during  the  computations  and  is  geared  to  deal  reliably  with  a  large 
class  of  problems.  Because  this  problem  is  so  simple,  the  100%  surcharge  is 
irrelevant . 

The  adaptively  constructed  grids  in  one  of  the  runs  are  pictured  in  Figure 
7.  The  input  and  results  for  this  run  are  summarized  in  the  row  marked  with 
an  asterisk  (*)  in  Table  II.  The  construction  of  the  grid  at  t  »  .6  was 
described  in  Section  7  and  pictured  in  Figure  5.  We  see  that  the  nodes  lead 
the  spreading  of  {uxx»vxx^  anc*  l°ca-*-  errors.  The  selection  of  these 

grids  was  "guided"  with  the  3-level  macro  grid  A(8,16,32).  The  parameter 
EPSDN,  allowed  in  all  examples  to  vary  in  [.6  EPS,  9  EPS],  was  chosen  here  to 
be  no  less  than  .7  EPS. 

Another  set  of  grids  for  this  problem  are  pictured  in  Figure  8.  These 
grids  were  constructed  by  the  adaptive  FEMOL  with  the  parameter  EPSDN  fixed 
as  .9  EPS  -  i.e.  the  method's  freedom  to  determine  the  grids'  intensities 
(cf.  Section  5)  was  restricted,  as  the  relative  | | | • | | | -error  was  lowered  to 
just  under  EPS  in  each  regridding.  The  number  of  nodes  used  in  this  run  was 
smaller  than  in  that  corresponding  to  Figure  7,  but  the  cost  as  measured  in 
elapsed  CPU  time  for  this  latter  run  was  60%  higher,  while  the  maximum 


re lative 


•  -error  was  no  lower. 
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TABLE  I. 

Results  for  Example  1  (NONADAPTIVE) :  REL  ERR  -  |  |  |e|  |  |  /  |  |  fu|  |  | 

(  u  ■  (u,v)  )  and  9  -  E  /  |||e|||  conputed  at  tp  -  p  ;  p  •  1 , 10  . 


NO.  ELTS. 

ODE  TOL 

REL 

ERR 

|e-i| 

CPO 

(unif . } 

stx 

ctan  max 

*p 

s 

lp  ep 

32 

tn 

1 

o 

—9 

• 

.19 

.20 

.016  .030 

1 

64 

t  /  2 

.09 

.10 

.003  .005 

2 

128 

T  /  4 

.05 

.05 

- 

6 

236 

T  /  8 

.02 

.02 

21 

TABLE  II. 

Reaulta  for  Example  1  (ADAPTIVE) :  REL  ERR  6  0  computed  aa  in  TABLE  I, 
32  unlf.  ELTS.  initially,  3  -  level  macro  grid  4(8,16,32)  uaed, 

REL  ERR  estimated  &  regridding  permitted  at  30  times  in  (0,10]. 
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NO. 
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ODE 
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|e-i 

U 

CPU 
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lP 
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2 

40 

44 
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.025 
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- 
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Example  2.  FltzHugh-Nagumo  Equations  (cf.  [18],  [21]) 

The  version  of  these  equations  considered  here  provides  a  conceptual  model 
of  ionic  current  flow  across  a  semi- inf inite  nerve  membrane.  An  electro¬ 
chemical  potential  u  and  "recovery"  variable  v  satisfy 


and 


where 


and  here 


u  ■  u  +  f(u)  -  v  ^ 
t  XX  I 


V  » 

t 


b(u  -  cv) 


t  >  0 

x  €  (0,120) 


(34) 


ux(t,0)  =*  -1/2 

ux(t,120)  -  0 


u  «  v  »  0  t  *  0  (36) 

x  €  (2,120) 

f(u)  *  u(u-a)(l-u)  (37) 

a  *  .139 

b  -  .008 

c  ■  ’  2 . 54 
I  -  .45 


(38) 


t  >  0 


(35) 


I  represents  the  magnitude  of  a  constant  current  applied  at  the  left  end 
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of  Che  nerve  and  b  is  Che  reciprocal  of  Che  Cime  scale  associaced  wich  Che 
nerve's  recovery.  Numerical  sCudies  in  Micchell  and  Manoranjan  [18],  Rinzel 
[21]  and  elsewhere  indicaCe  Che  sensicivicy  of  Che  solucion  behavior  Co 
changes  in  Che  paramecers  a,  b,  c,  and  I.  Preliminary  sCudies  here 
suggesc  chaC  wich  Che  values  in  (38),  repeCiCive  pulses  (in  u  and  v) 
Craveling  aC  speed  *  .4  are  generaced  wich  firing  frequency  «  .77  x  10*^. 
Firing  frequency  is  defined  as  Che  reciprocal  of  Che  Cemporal  period  for  a 
solucion  appearing  Co  be  a  craveling  wave  for  large  x  and  c. 

The  evolucion  of  u  *  {u,v}  is  shown  in  Figure  9,  where  values  were 
compuCed  in  one  of  Che  mosc  accuraCe  nonadapcive  runs.  Ic  is  noced  ChaC 
because  of  Che  diffusionless  form  of  Che  second  of  Eqs.  (34),  only  ux  -  Ux 
is  expliciCly  concrolled  in  Che  adapCive  FEMOL,  and  Che  accuracy  of 
approximaCing  v  encers  naCurally  Chrough  ics  role  in  Che  residual  of  Che 
firsC  of  Eqs.  (34). 

The  resulCs  for  chis  problem  are  summarized  in  Tables  III  and  IV.  The 
relacive  I t | *  I  I | —error  and  0  were  compuced  wich  Che  "exacc"  solucion  as 
obCained  via  960  uniform  elemencs  and  very  accuraCe  ODE  Cime  inCegraCion  in 
Che  nonadapcive  FEMOL.  EscimaCed  values  were  used  for  Che  relacive  errors  in 
Che  mosC  accuraCe  encries  of  chese  Cables,  as  was  done  in  Example  1. 

The  qualiCy  of  E  is  seen  in  Che  | 8—1 |  columns  of  Tables  III  and  IV. 

AC  Cime  C  *  80,  when  a  new  pulse  is  rapidly  developing  ac  Che  lefc  spacial 
boundary,  tie  error  was  underesCimaCed  by  43 7.  when  a  large  (.2)  relacive 
error  was  requesCed  in  Che  adapCive  FEMOL.  E  esCimaCed  Che  error  much  more 
accuraCely  ac  Che  ocher  cargec  Cimes,  however,  and  Che  qualiCy  of  E 
evidenCly  improves  as  Che  relacive  j | | • | I | —error  ♦  0,  in  boCh  Che  adapCive 
and  nonadapcive  cases. 

By  comparing  Che  EPS  and  max  REL  ERR  encries  of  Table  IV,  two  failures  Co 
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achieve  error  control  according  to  the  goal  (11)  are  seen,  the  failure  for  EPS 
*  ,1  being  "lower  level"  than  that  for  EPS  *  .2.  These  failures  could  have 
been  eliminated  by  "tuning"  the  algorithm — i.e.g.  allowing  regridding  to  occur 
before  E/|||u|||  >  .95  EPS  -  or  by  requesting  more  accuracy  than  is  actually 
needed .  But  these  failures  show  how  dependent  goal  achievement  is  on  the 
performance  of  the  error  estimator  in  an  adaptive  method  such  as  the  present. 

By  comparing  the  NO.  ELT.  and  CPU  entries  of  the  last  three  rows  of  Tables 
III  and  IV,  one  sees  that  the  mean  number  of  elements  used  in  the  adaptive 
FEMOL  was  about  Vo  of  the  corresponding  number  of  uniform  elements,  but  that 
the  same  fraction  of  CPU  time  was  not  saved  until  the  l | I • | l | —error  was 
about  2-3Z.  This  is  primarily  due  to  the  fact  that  the  traveling  fronts  in 
this  problem  are  not  very  localized.  As  for  Example  1,  we  emphasize  that  the 
CPU  times  for  the  nonadaptive  runs  do  not  reflect  the  true  cost  of  achieving 
EPS  relative  | J j • | | | -accuracy,  whereas  those  for  the  adaptive  runs  do. 

Note  that  the  number  of  regriddings  remained  relatively  constant  as  EPS 
and  TOL  ♦  0  and  that  the  mean  number  of  finite  elements  appears  to  have 
linearly  increased  with  decreasing  EPS.  These  observations  are  made  with 
almost  all  of  the  experiments  of  this  section.  They  suggest  two  things: 

(1)  The  regridding  strategy  is  being  carried  out  exactly  as  planned. 

Grid  shapes  are  chosen  (somewhat  stablely)  in  adapting  to  the 
changing  nature  of  the  solution,  while  the  numbers  of  nodes  (or  grid 
intensities)  are  determined  by  the  requested  accuracy  EPS. 

(2)  Efficiency  demands  that  pattern  recognition  notions  be  used  to 
predict  grid  shapes,  since  the  number  of  regriddings  will  not 
increase  unboundedly  as  EPS  +  TOL  ♦  0.  The  workhorse  of  the  present 
method  is  the  implicit  ODE  solver  LSODI.  The  time  stepsizes  used  by 
LSODI  for  each  space  grid  decrease  as  the  error  tolerances  decrease. 
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but  the  relative  distributions  (or  shapes)  of  the  stepsize  sequences 
change  very  little.  These  distributions  depend  primarily  on  the 
shapes  of  the  space  grids  and  higher  solution  derivatives,  and  they 
govern  the  most  appropriate  length  of  time  to  retain  a  grid. 

The  grids  constructed  in  one  of  the  runs  for  this  problem  are  shown  in 
Figure  10.  The  input  and  results  for  this  run  are  summarized  in  the  row 
marked  with  an  asterisk  (*)  in  Table  IV.  Here,  regridding  occurred  at  20  of 
the  100  times  at  which  it  was  allowed.  Grid  shape  selection  was  guided  with 
the  3- level  macro  grid  4(10,20,40),  and  the  parameter  EPSDN  determining  grid 
intensity  was  chosen  to  be  no  less  than  .8  EPS. 
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TABLE  III. 

R*«ult*  for  Exaapla  2  (NONADAPTIVB) :  EEL  ERR  -  (((«[[(  /  | | |u | | | 


(u,v)  )  and  •  • 

E  /  II 

.Mil 

coaputad 

at 

tp  ■  2  Op  ; 

P  "  1,4 , 

NO.  ELTS. 

ODE 

TOL 

REL 

ERR 

Ie-I| 

CPU 

(unlf . ) 

■•an 

■an 

■aan  max 
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P  P 
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8 
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.OA  .08 

25 
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A 

.05 
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73 
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8 
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TABLE  IV. 

Raaulta  for  Exaapla  2  (ADAPTIVE):  REL  ERR  &  8  coaputad  aa  In  TABLE  III, 
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Example  3 .  Two  Counter-Traveling  Waves 

This  model  problem  was  created  in  order  to  examine  the  adaptive  FEMOL's 
performance  when  two  wave  forms,  having  different  front  widths  and  directions 
and  speeds  of  travel,  collide  and  pass  through  one  another.  The  solution  of 
this  problem  is 


U  -  ud)  +  u<2> 


(39) 


where 


u(1)(t,x)  -  .25(1  +  tanh( 100(x-10t))) 

u(2)(t,x)  =  .25(1  +  tanh(80( l-x-30t) )) 


(40) 


u<*>  moves  to  the  right  at  speed  10  and  u^2\  whose  front  width  is  25% 
larger  than  that  of  moves  to  the  left  at  speed  30.  At  time  t  =  .025 

the  waves  collide,  almost  extinguishing  one  another.  The  PDE  numerically 
solved  was  the  scalar  heat  equation 


ut  "  “xx  +  f  *  t  >  0,  x  €  (0,1)  (41) 

where  f,  the  initial  data  at  t  *  0,  and  Dirichlet  boundary  conditions  were 
set  according  to  (39)  and  (40). 

The  results  for  this  problem  are  summarized  in  Tables  V  and  VI.  The  high 
quality  and  rapid  convergence  of  the  error  estimator  E  are  suggested  by  the 
| 0-1 |  columns  of  both  tables.  One  notices  that  some  of  this  quality  is 
sacrificed  when  the  grids  are  nonuniform  and  allowed  to  change.  The 
performance  of  E  here  is  typical  of  that  we  have  seen  in  linear  and  most 
mildly  nonlinear  problems,  in  which  the  nonlinearities  are  polynomial  or 


rational  in  form 
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Error  control  was  achieved  In  all  adaptive  runs,  with  the  exception  of  the 
low  level  failure  that  occurred  with  EPS  *  .05.  One  sees  by  comparing  CPU 
times  in  Tables  V  and  VI  that  the  cost  in  achieving  this  was  only  less  than 
that  for  the  uniform,  nonadaptive  case  when  the  error  was  about  2-3%.  This 
cost  comparison  of  course  excludes  the  overhead  associated  with  determining 
the  number  of  uniform  elements  needed  to  obtain  a  given  relative 
| | ) • | | | -accuracy  EPS. 

The  grids  constructed  in  one  of  the  adaptive  runs  are  pictured  in  Figure 
11,  along  with  the  locations  of  the  front  centers.  The  input  and  results  for 
this  run  are  summarized  in  the  row  marked  with  an  asterisk  (*)  in  Table  VI. 

One  sees  that  for  0  <  t  <  .025,  the  regridding  frequency  was  determined 
automatically  by  the  faster  left-moving  front,  and  that  nodes  were  placed 
ahead  of  it  in  an  effective  manner.  For  t  >  .035  the  same  was  true  of  the 


right-moving  front 
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TABLE  V. 

Result*  for  Example  3  (NONADAPTIVE) :  REL  ERR  -  |||*|||  /  |||u||| 
&  0  -  E  /  111*111  computed  at  tp  -  .01  p  ;  p  •  1,7  . 
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TABLE  TI. 

Results  for  Example  3  (ADAPTIVE) :  REL  ERR  &  6  computed  as  In  TABLE  V, 
80  unif.  ELTS.  Initially,  3  -  level  macro  grid  6(10, 20,40)  used, 

REL  ERR  estimated  &  regrlddlng  permitted  at  140  times  In  (0..07J  . 
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FIN.  aTS. 


Figure  11.  Adaptively  constructed  grids  for  Example  3 


and  wave  centers  x 


10  t 


and  x 


1 


30  t  . 
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Example  4.  Dwyer-Sanders  Model  Flame 

This  model  has  been  proposed  to  simulate  many  of  the  basic  features  of 
flame  propagation.  It  was  used  as  a  test  problem  for  the  Moving  Finite 
Element  method  in  Gelinas,  et  at.  [12].  The  equations  governing  a  single 
species  mass  density  p  and  temperature  T  are 


and 


where 


pt  "  pxx  -  f(T)p 
Tt  “  Txx  +  f<T>P 


t  >  0 

[  x  6  (0,1) 


(42) 


Px(t,0) 

Tx(t,0) 


Px(t,l)  -  0 


0,  T(t , 1)  *  g(t) 


y  t  >  o 


(43) 


P  -  1,  T  -  .2 


t  -  0,  (44) 

x  e  (0,1) 


and 


f(T)  -  3.52  x  io6  exp(-4/T) 


(45) 


g(t) 


.2  +  t/(2  x  10-4)  ; 

.1.2  ; 


t  <  2  x  10-4, 
t  >  2  x  10-4. 


(46) 


The  flame  is  ignited  at  (t,x)  -  (0,1)  and  propagates  from  right  to  left 
at  a  relatively  high  rate  of  speed.  In  the  numerical  studies  here,  p  and 
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T  were  computed  very  accurately  for  0  <  t  <  2  x  10  \  and  the  values  at  t 

■  2  *  10*4  were  used  as  initial  conditions  for  all  experimental  runs,  in 

which  integrations  were  carried  out  until  time  t  «■  .006. 

The  steep  temperature  profiles  in  the  moving  temperature-mass  density 
front  which  were  computed  in  one  of  the  most  accurate  runs  are  pictured  in 

Figure  12.  The  front  propagates  at  approximately  constant  speed  for  .003  <  t 

<  .006.  This  problem  is  much  more  difficult  than  the  others  considered  thus 
far,  due  to  the  nature  of  the  nonlinear  kinetics.  Changing  grids  to  control 
errors  and  accurately  resolve  flame  structure  is  more  of  a  necessity  here  than 
a  convenience. 

The  results  of  applying  the  adaptive  FEMOL  to  this  problem  with  various 
requested  relative  accuracies  EPS  are  summarized  in  Table  VII.  The 
experimental  procedure  differed  from  that  used  in  Examples  1-3  in  two  ways. 
First,  it  was  not  feasible  to  obtain  an  almost  "exact"  solution  for  comparison 
by  using  uniform,  nonadaptive  finite  elements.  The  first  three  entries  in  the 
REL  ERR  and  | 0—1 |  columns  of  Table  VII  were  therefore  computed  by  using  the 
approximate  solution  obtained  with  the  adaptive  FEMOL  and  the  parameters  in 
row  5  as  "exact".  The  estimated  accuracy  of  this  approximate  solution  was 
better  at  many  of  the  12  target  times  than  is  indicated  by  the  estimated  mean 
and  max  REL  ERR  entries  in  row  5.  Second,  the  same  small  input  ODE  error 
tolerance  T0L  was  used  in  each  run.  This  was  done  because  the  dependence  of 
the  time  discretization  error  on  T0L  was  observed  to  be  less  smooth  in  this 
problem  than  in  the  others . 

The  approximate  wave  speed  for  .003  <  t  <  .006  in  each  of  the  runs  is 
listed  in  the  last  column  of  Table  VII.  The  speeds  were  computed  by 
monitoring  the  spatial  locations  of  the  T  ■  .5,  .75,  and  1.0  values  on 
the  moving  temperature  front  at  7  evenly  spaced  times.  Those  familiar  with 


this  model  will  recognize  the  accuracies  of  these  speeds,  and  rapid 
convergence  as  EPS  ♦  0  is  suggested  by  the  entries. 

With  the  larger  input  EPS  values,  the  ability  of  E  to  estimate  the 
111*11 |— error  decays  significantly  as  time  increases  past  .003.  This  can  be 
seen  by  comparing  each  of  the  mean  and  max  | 9—1 |  values  in  Table  VII  and 
would  be  expected,  given  the  relatively  inaccurate  speeds  of  propagation.  The 
quality  of  E  improves  as  REL  ERR  0,  however.  This  is  suggested  by  the 
third  pair  of  | 0—1 (  entries  and  the  mass  density  and  temperature  gradient 
profiles  computed  with  various  EPS.  The  temperature  gradients  computed  in 
many  of  the  runs  at  the  final  time  are  shown  in  Figure  13. 

The  ability  of  the  adaptive  scheme  to  control  errors  according  to  (11) 
depends  strongly  on  its  ability  to  estimate  the  errors.  High  level  failures 
to  achieve  (11)  are  seen  in  Table  VII  for  large  EPS  values  and  the  reason  for 
these  failures  is  evident  from  the  similarity  of  the  REL  ERR  and  } 9—1 | 
entries.  We  believe  that  if  failures  occurred  for  the  two  smallest  EPS 
values,  however,  then  these  failures  were  at  a  much  lower  level  than  those 
observed  for  the  large  EPS  values . 

In  order  to  relate  the  cost  of  the  adaptive  FEMOL  to  that  incurred  with 
uniform,  nonadaptive-eleraents  for  this  problem,  1000  elements  and  the  same  TOL 
value  were  used  in  one  nonadaptive  run.  Over  1200  CPU  seconds  were  required 
to  complete  the  integration.  The  mean  wave  speed  (as  measured  before)  was 
142.19,  which  along  with  the  observed  mass  density  and  temperature  gradient 
profiles  indicates  that  the  same  accuracy  would  have  been  achieved  with  EPS  - 
.04  in  the  adaptive  FEMOL,  at  less  than  1/6  of  the  cost. 

The  automatically  constructed  grids  in  one  of  the  runs  are  shown  in  Figure 
14,  along  with  the  flame  center  locations  for  .003  <  t  <  .006.  The  input  and 
results  for  the  run  in  which  these  grids  were  constructed  are  summarized  in 
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the  row  marked  with  an  asterisk  (*)  in  Table  VII.  The  selection  of  the  space 
grids  was  guided  with  the  3-level  macro  grid  £(12,24,48).  The  construction 
of  the  grid  at  t  *  3,8  x  10-^  was  described  in  Section  7  and  pictured  in 
Figure  4.  Clearly,  nodes  are  distributed  in  a  way  which  anticipates  the 
moving  front. 


Figure  12. 


Temperature  T  at  various  times  in 
flame  model  in  Example  4. 
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TABLE  VII. 

Raaulta  for  Exaapla  A  (ADAPTIVE);  EEL  EXE  -  111*111  /  |||u||| 

(  u  *  { T,  p )  )  4  8  "  E  /  M  1*1  I  I  coaputad  *t  tp  ■  5p  x  10~*  ;  p  <•  1,12  . 

192  unif.  ELTS.  Initially,  3  -  laval  ncro  grid  6(12,24,48)  usad, 

ODE  TOL  •  2  x  10"6  In  all  rttna,  EEL  EXE  aatlaatad  4  ragrlddlng 
paralttad  at  290  tlaaa  In  (.0002  ,  .006]  . 


EPS 

MESH 

NO. 

ELTS. 

EEL 

EBB 

l»-1 

u 

CPD 

HAVE 

HODS 

Bean 

BAX 

aaan 

£P 

aax 

CP 

nun 

‘p 

aax 

*p 

SPEED 

e  -  .2 

12 

62 

103 

.55 

.90 

.64 

.86 

43 

143.6 

t  1 

.1 

*  E  /  2 

22 

111 

126 

.18 

.35 

.58 

.83 

80 

142.6 

+ 

.3 

t  /  4 

20 

170 

188 

.05 

.07 

.22 

.43 

133 

142.23 

+ 

.05 

c  /  8 

18 

289 

353 

.02 

.03 

- 

- 

275 

142.11 

.  02 

c  / 16 

20 

487 

724 

.01 

.01 

- 

- 

680 

142.08 

+ 

.01 
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We  conclude  by  describing  Che  results  of  an  experiment  designed  to  assess 
the  cost  effectiveness  of  decisions  made  in  adaptive  regridding.  The  most 
important  decisions  concern  grid  shape  and  intensity.  A  grid's  intensity  is 
controlled  by  selecting  the  parameter  EPSDN  6  [.6  EPS,  9  EPS]  (cf.  Sections 
5,  6).  A  grid's  shape  is  controlled  by  choosing  1  of  3  levels  as  being 
"active"  on  each  of  the  J  level  1  macro  elements  contained  in  A(J,2J,4J) 
(cf.  Sections  5,  7). 

In  each  of  the  four  examples  of  this  section  we  restricted  the  algorithm's 
freedom  to  choose  these  parameters  and  compared  costs  with  that  for  an 
"unrestricted"  run.  The  unrestricted  run  in  each  example  is  that 
corresponding  to  the  *  table  entries  and  that  which  produced  the  grids 
pictured.  Eleven  restricted  runs  were  made  for  each  problem:  EPSDN  was  fixed 
as  .7  EPS,  .8  EPS  and  .9  EPS,  the  "active"  levels  in  all  level  1  macro 
elements  were  fixed  as  the  largest  and  the  smallest  possible ,  and  all 
combinations  of  these  EPSDN  and  "active"  level  restrictions  were  imposed. 

The  CPU  times  used  in  the  44  restricted  parameter  runs  were  rounded  to  the 
nearest  multiple  of  5X  of  the  corresponding  unrestricted  run  CPU  times.  These 
results  are  summarized  by  the  histogram  in  Figure  15.  While  some  of  the 
restricted  runs  were  slightly  less  expensive,  we  see  that  the  decisions  made 
by  the  algorithm  are  lowering  cost  the  vast  majority  of  the  time. 


NUM.  OF 
RUNS 


Figure  15 
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1  CPU  TIME  OVER  CORRESPONDING 
UNRESTRICTED  PARAMETER  RUNS 


Added  cost  in  restricting  the  adaptive 
algorithm's  freedom  in  the  four  examples 


of  this  section. 
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9.  SUMMARY 

Using  piecewise  linear  finite  elements  in  one  space  dimension  and  implicit 
integration  formulas  in  time,  an  adaptive  method  of  lines  has  been  developed 
for  a  class  of  nonlinear  parabolic  partial  differential  equations.  This  class 
includes  many  equations  which  are  used  to  model  reaction-diffusion  phenomena. 

The  primary  goal  of  this  adaptive  method  is  to  keep  a  particular  norm  of 
the  space  discretization  error  less  than  a  user-specified  tolerance  EPS  at  all 
times.  Error  control  is  obtained  by  adding  and  deleting  space  nodes  when  a 
computed  estimate  of  the  error  has  nearly  exceeded  EPS.  This  process  is  made 
efficient,  without  raising  the  risk  of  an  error  control  failure,  by  monitoring 
ODE  time  stepsize  information  and  utilizing  multi-grid  pattern  recognition 
notions  to  predict  appropriate  grid  intensities  and  shapes. 

Experiments  presented  here  have  shown  that  the  computed  error  estimates 
are  accurate,  that  the  procedure  reliably  controls  errors,  and  that  the 
strategy  to  keep  costs  low  is  successful.  These  experiments  were  conducted 
with  the  program  FEM0L1,  which  uses  the  LSODI  package  to  solve  the  ODE  initial 
value  problem  resulting  from  each  space  grid.  Information  concerning  FEM0L1 
can  be  obtained  by  contacting  the  first  author. 

The  process  of  collecting  information  about  a  solution  through  local 
estimates,  assimilating  this  feedback  at  a  more  global  level,  and  utilizing 
features  extracted  from  the  global  representations  was  shown  here  to  be  an 
effective  and  efficient  way  to  automatically  construct  grids.  Many  adaptive 
schemes  incorporate  this  process  in  one  form  or  another,  but  the  question  of 
how  this  should  best  be  done  has  not  been  adequately  addressed. 
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